******************************************************************************
*										
*                           Boese & Eberhardt							
*            
*                          Facets of Democracy
*													
******************************************************************************
*                           Descriptives 2024
******************************************************************************
*                           Created: 6th May 2021					
*                      Last Changed: 17th July 2024
******************************************************************************
*                              Markus Eberhardt						
******************************************************************************
*                 School of Economics, University of Nottingham,			
*                  University Park, Nottingham NG7 2RD, England			
*                     email: markus.eberhardt@nottingham.ac.uk				
*                   web: https://sites.google.com/site/medevecon/			
******************************************************************************

clear all
set matsize 11000

global path "D:/Dropbox/Facets of Democracy"
global rawpath "D:/Dropbox"
global otherpath "/Users/lezme/Dropbox"
global path "/Users/lezme/Dropbox/Facets of Democracy"
global rawpath "/Users/lezme/Dropbox"


global data "$path/Data"
global otherdata "$path/Data"
global do "$path/Do_files"

global texfolder "$rawpath/Apps/Overleaf/Facets of Democracy"


****
**** Get merged dataset and prepare variables
****

use "$data/madd_ert_vdem_2021.dta", clear
tsset nwbcode year

* Maddison: per capita GDP (in logs)*100 and population growth rate
gen y = ln(gdppc)*100
gen lpop = ln(pop)
gen dlpop = 100*d.lpop

* DOTS data: export/trade
gen ex = ex_trade*100

xtreg y dlpop ex v2x_libdem if year>=1959, fe
* n=8,303, N=157, avg T=52.9 (min 12, max 70); 1949-2018

drop csum 
gen c=1 if e(sample)
by nwbcode: egen csum=sum(c) if c==1
gen sample= (e(sample))

* ROW: democracy dummy
gen dem_row2 = (v2x_regime>=2) & !missing(v2x_regime)
gen dem_row2amb = (v2x_regime_amb>=6) & !missing(v2x_regime)
label var dem_row2 "ROW - Regime Change 1 to 2"

* PolityIV: democracy dummy
gen dem_p4_1 = (polity2>=1) & !missing(polity2)
gen dem_p4_6 = (polity2>=6) & !missing(polity2)

* V-Dem: democracy dummies from high and mid-level indicators
local z "libdem polyarchy liberal"
foreach k of local z{
	sum v2x_`k' if sample==1
	scalar mean`k' = r(mean) 
	scalar sd`k' = r(sd) 
	gen st_`k' = (v2x_`k'-mean`k')/sd`k' if sample==1
}

gen dem_ld = (st_libdem>0)
gen dem_poly = (st_polyarchy>0)
gen dem_lib = (st_liberal>0)

* V-Dem: democracy dummies from high and mid-level indicators
local z "_freexp_altinf _frassoc_thick el_frefair cl_rol _jucon lg_legcon"
foreach k of local z{
	sum v2x`k' if sample==1
	scalar mean`k' = r(mean) 
	scalar sd`k' = r(sd) 
	gen st`k' = (v2x`k'-mean`k')/sd`k' if sample==1
}

* Positive values
local z "0"
foreach k of local z{
	local norm_`k' = `k'/100
	gen dem_poly_`k' = (st_polyarchy>=`norm_`k'') & !missing(v2x_polyarchy)
	gen dem_lib_`k' = (st_liberal>=`norm_`k'') & !missing(v2x_liberal)
	gen dem_exp_`k' = (st_freexp_altinf>=`norm_`k'') & !missing(v2x_freexp_altinf)
	label var dem_exp_`k' "Freedom of Expression Threshold > mean + .`k' sd"
	gen dem_assoc_`k' = (st_frassoc_thick>=`norm_`k'') & !missing(v2x_frassoc_thick)
	label var dem_assoc_`k' "Freedom of Association Threshold > mean + .`k' sd"
	gen dem_fairel_`k' = (stel_frefair>=`norm_`k'') & !missing(v2xel_frefair) 
	label var dem_fairel_`k' "Clean Election Threshold Threshold > mean + .`k' sd"
	gen dem_rolaw_`k' = (stcl_rol>=`norm_`k'') & !missing(v2xcl_rol)
	label var dem_rolaw_`k' "Rule of Law Threshold > mean + .`k' sd"
	gen dem_jucon_`k' = (st_jucon>=`norm_`k'') & !missing(v2x_jucon)
	label var dem_jucon_`k' "Jud. Constraints Threshold > mean + .`k' sd"
	gen dem_legcon_`k' = (stlg_legcon>=`norm_`k'') & !missing(v2xlg_legcon) 
	label var dem_legcon_`k' "Leg. Constraints Threshold > mean + .`k' sd"
}

* Cheibub et al  
rename democracy dem_cheibub 
label var dem_cheibub "Cheibub et al Democracy indicator (1950-2008)"
label var dem_boix "Boix et al Democracy indicator (1950-2014)"
label var dem_ANRR "ANRR Democracy indicator (1960-2010)"

save "$data/madd_ert_vdem_2024_helper2.dta", replace

***
*** 'Birth' of a country (as far as the sample is concerned)
***

gen data = 1 if !missing(v2x_libdem) & sample==1
sort nwbcode year
by nwbcode: gen sumdata = sum(data) 

********************************************************************************
*** INCLUDED IN THE DRAFT: When did countries appear in the dataset?
********************************************************************************
twoway ///
	(hist year if sumdata==1, w(10) fraction fcolor(midblue*.5%50) lw(none)) ///
	(hist year if sumdata==1, w(5) fraction fcolor(midblue%50) lw(none)) ///
	(hist year if sumdata==1, w(1) fraction  lcolor(pink*1.25) fcolor(pink*1.25%50) lw(.2)) ///
	(hist year if sumdata==1, w(10) fraction fcolor(none) lw(none) yaxis(2)) ///
	(hist year if sumdata==1, w(5) fraction fcolor(none) lw(none) yaxis(2)) ///
	(hist year if sumdata==1, w(1) fraction  lcolor(none) fcolor(none) lw(none) yaxis(2)) ///
, scheme(s1color) graphregion(color(white)) yscale(range(-0.005 .7001))  ymtick(##10) ///
	ysize(4) xsize(8) xlabel(1959(5)2009, labsize(small)) xscale(range(1958.01 2009.99)) xmtick(##5)  ///
	legend(label(1 "Decadal") label(2 "Quinquennial")  label(3 "Annual") rows(1) order(1 2 3)) ///
	ylabel(0(.10).7, labsize(small) angle(0) axis(2)) ytitle("Share of Countries", axis(2)) ///
	ylabel(0(.10).7, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash) angle(0)) ///
	ytitle("Share of Countries") xtitle("Year country first appears in the dataset (1949-2018)") ///
		yscale(range(-0.005 .7001) axis(2))  ymtick(##10, axis(2))
graph export "$texfolder/year_appear2024.eps", as(eps) replace
********************************************************************************
*** INCLUDED IN THE DRAFT
********************************************************************************

xtlogit dem_poly y if sample==1, fe
gen data2 = 1 if e(sample)
sort nwbcode year
by nwbcode: gen sumdata2 = sum(data2) 

********************************************************************************
*** INCLUDED IN THE DRAFT: When did countries appear in the dataset? (POLY)
********************************************************************************
twoway ///
	(hist year if sumdata2==1, w(10) fraction fcolor(midblue*.5%50) lw(none)) ///
	(hist year if sumdata2==1, w(5) fraction fcolor(midblue%50) lw(none)) ///
	(hist year if sumdata2==1, w(1) fraction  lcolor(pink*1.25) fcolor(pink*1.25%50) lw(.2)) ///
	(hist year if sumdata2==1, w(10) fraction fcolor(none) lw(none) yaxis(2)) ///
	(hist year if sumdata2==1, w(5) fraction fcolor(none) lw(none) yaxis(2)) ///
	(hist year if sumdata2==1, w(1) fraction  lcolor(none) fcolor(none) lw(none) yaxis(2)) ///
, scheme(s1color) graphregion(color(white)) yscale(range(-0.005 .7001)) ymtick(##10) ///
	ysize(4) xsize(8) xlabel(1959(5)2009, labsize(small)) xscale(range(1958.01 2009.99)) xmtick(##5)  ///
	legend(label(1 "Decadal") label(2 "Quinquennial")  label(3 "Annual") rows(1) order(1 2 3)) ///
	ylabel(0(.10).7, labsize(small) angle(0) axis(2)) ytitle("Share of Countries", axis(2)) ///
	ylabel(0(.10).7, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash) angle(0)) ///
	ytitle("Share of Countries") xtitle("Year country first appears in the dataset (1949-2018)") ///
	yscale(range(-0.005 .7001) axis(2))  ymtick(##10, axis(2))
graph export "$texfolder/year_appear2024_poly.eps", as(eps) replace
********************************************************************************
*** INCLUDED IN THE DRAFT
********************************************************************************


***
*** Histograms
***

use "$data/madd_ert_vdem_2024_helper2.dta", clear

gen demevent_ld = 1 if ((dem_ld==1 & l.dem_ld==0)) & sample==1
gen demevent_poly = 1 if ((dem_poly==1 & l.dem_poly==0)) & sample==1
gen demevent_lib = 1 if ((dem_lib==1 & l.dem_lib==0)) & sample==1

gen demevent_exp = 1 if ((dem_exp_0==1 & l.dem_exp_0==0)) & sample==1
gen demevent_assoc = 1 if ((dem_assoc_0==1 & l.dem_assoc_0==0)) & sample==1
gen demevent_fairel = 1 if ((dem_fairel_0==1 & l.dem_fairel_0==0)) & sample==1
gen demevent_rolaw = 1 if ((dem_rolaw_0==1 & l.dem_rolaw_0==0)) & sample==1
gen demevent_jucon = 1 if ((dem_jucon_0==1 & l.dem_jucon_0==0)) & sample==1
gen demevent_legcon = 1 if ((dem_legcon_0==1 & l.dem_legcon_0==0)) & sample==1


twoway ///
	(hist year if demevent_ld==1 & sample==1, discrete frequency lcolor(gold) fcolor(orange*1.25%50) lwidth(.2) barwidth(1)) ///
	(hist year if demevent_poly==1 & sample==1, discrete frequency lcolor(gs6) fcolor(pink*1.25%50) lwidth(.2) barwidth(.75)) ///
	(hist year if demevent_lib==1 & sample==1, discrete frequency lcolor(gs6) fcolor(midblue*1.25%50) lwidth(.2) barwidth(.5)) ///
, legend(order(1 2 3) rows(1) ///
	label(1 "[113 events]") label(2 " [61 events]") label(3 " [61 events]")) /// 
	scheme(s2mono) graphregion(color(white)) yscale(range(0 8.01)) ///
	ytitle("Count (unbalanced panel)") ylabel(0(1)8, grid glc(black) glw(.2) glp(shortdash)) ///
	xlabel(1950(5)2015) ymtick(##5)	tmtick(##5) ysize(4) xsize(9) xscale(range(1949 2019))
	
tab demevent_ld	

***
*** Scatter plots (Single difference estimates)
***

gen lngdppc = ln(gdppc)

*** POLYARCHY
xtlogit dem_poly_0 y dlpop ex if sample==1, fe
sort nwbcode year
by nwbcode: egen dem_sum=sum(sample) if e(sample) & dem_poly_0==1
by nwbcode: egen end_sample=max(year) if e(sample) & dem_poly_0==1
by nwbcode: egen start_sample=min(year) if e(sample) & dem_poly_0==1
gen yT = ln(gdppc) if year==end_sample
by nwbcode: egen yT_ = mean(yT)
gen y0 = ln(gdppc) if year==start_sample
by nwbcode: egen y0_ = mean(y0)
gen changey = 100*(d.lngdppc) if e(sample) & dem_assoc_0==0
by nwbcode: egen change__y = mean(changey)
gen change_y = ((100*(yT_ - y0_)/dem_sum)/(change__y))*dem_sum

sort nwbcode year
by nwbcode: egen sum_poly = sum(demevent_poly) if sample==1
tab sum_poly if  year==end_sample

twoway ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y>-100 & change_y<150 & sum_poly==1, ///
		msymbol(Oh) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mlw(medium) mcolor(gs5)) ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y>-100 & change_y<150 & sum_poly==2, ///
		msymbol(Oh) msize(medlarge)  mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mlw(thick) mcolor(gs5)) ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y>-100 & change_y<150 & sum_poly>=3, ///
		msymbol(O) msize(medlarge)  mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mcolor(gs5)) ///
	(fpfit change_y dem_sum if e(sample) & year==end_sample, ///
		lcolor(pink*1.25) lw(.8) lpattern(solid)) ///
, xsize(11) ysize(6) scheme(s2mono)  graphregion(color(white) margin(small)) plotregion(margin(small)) ///
	xlabel(0(5)70, labsize(small)) yscale(range(-51 151)) ///
	ylabel(-50(50)150, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash) angle(0)) ///
	xtitle("Years Spent in Regime (Length of Treatment)", size(small)) ///
	ytitle("Total Change in Real GDP pc during Regime (in %)", size(small))  ///
	yline(0, lpattern(shortdash) lw(.4) lcolor(black)) xmtick(##5) ///
	text(125 0 "{bf:Polyarchy} (mean threshold)", size(medsmall) place(e)) ///
	legend(order(- "Number of Regime Changes" 1 2 3 4) label(1 "1") label(2 "2")  label(3 "3-4") ///
		label(4 "Fractional Polynomial") rows(1) size(small) ) 
graph export "$texfolder/Drilling_scatter_poly.eps", as(eps) replace

drop dem_sum end_sample yT start_sample y0 yT_ y0_ change_y changey change__y


*** LIBERAL COMPONENT
xtlogit dem_lib_0 y dlpop ex if sample==1, fe
sort nwbcode year
by nwbcode: egen dem_sum=sum(sample) if e(sample) & dem_lib_0==1
by nwbcode: egen end_sample=max(year) if e(sample) & dem_lib_0==1
by nwbcode: egen start_sample=min(year) if e(sample) & dem_lib_0==1
gen yT = ln(gdppc) if year==end_sample
by nwbcode: egen yT_ = mean(yT)
gen y0 = ln(gdppc) if year==start_sample
by nwbcode: egen y0_ = mean(y0)
gen changey = 100*(d.lngdppc) if e(sample) & dem_assoc_0==0
by nwbcode: egen change__y = mean(changey)
gen change_y = ((100*(yT_ - y0_)/dem_sum)/(change__y))*dem_sum

sort nwbcode year
by nwbcode: egen sum_lib = sum(demevent_lib) if sample==1
tab sum_lib if  year==end_sample

twoway ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y>-100 & change_y<150 & sum_lib==1, ///
		msymbol(Oh) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mlw(medium) mcolor(gs5)) ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y>-100 & change_y<150 & sum_lib==2, ///
		msymbol(Oh) msize(medlarge)  mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mlw(thick) mcolor(gs5)) ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y>-100 & change_y<150 & sum_lib>=3, ///
		msymbol(O) msize(medlarge)  mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mcolor(gs5)) ///
	(fpfit change_y dem_sum if e(sample) & year==end_sample, ///
		lcolor(midblue*1.25) lw(.8) lpattern(solid)) ///
, xsize(11) ysize(6) scheme(s2mono)  graphregion(color(white) margin(small)) plotregion(margin(small)) ///
	xlabel(0(5)70, labsize(small)) yscale(range(-51 151)) ///
	ylabel(-50(50)150, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash) angle(0)) ///
	xtitle("Years Spent in Regime (Length of Treatment)", size(small)) ///
	ytitle("Total Change in Real GDP pc during Regime (in %)", size(small))  ///
	yline(0, lpattern(shortdash) lw(.4) lcolor(black)) xmtick(##5) ///
	text(125 0 "{bf:Liberal Component} (mean threshold)", size(medsmall) place(e)) ///
	legend(order(- "Number of Regime Changes" 1 2 3 4) label(1 "1") label(2 "2")  label(3 "3-6") ///
		label(4 "Fractional Polynomial") rows(1) size(small) ) 
graph export "$texfolder/Drilling_scatter_lib.eps", as(eps) replace

drop dem_sum end_sample yT start_sample y0 yT_ y0_ change_y change__y changey


*** FREEDOM OF EXPRESSION
xtlogit dem_exp_0 y dlpop ex if sample==1, fe
sort nwbcode year
by nwbcode: egen dem_sum=sum(sample) if e(sample) & dem_exp_0==1
by nwbcode: egen end_sample=max(year) if e(sample) & dem_exp_0==1
by nwbcode: egen start_sample=min(year) if e(sample) & dem_exp_0==1
gen yT = ln(gdppc) if year==end_sample
by nwbcode: egen yT_ = mean(yT)
gen y0 = ln(gdppc) if year==start_sample
by nwbcode: egen y0_ = mean(y0)
gen changey = 100*(d.lngdppc) if e(sample) & dem_assoc_0==0
by nwbcode: egen change__y = mean(changey)
gen change_y = ((100*(yT_ - y0_)/dem_sum)/(change__y))*dem_sum

sort nwbcode year
by nwbcode: egen sum_exp = sum(demevent_exp) if sample==1
tab sum_exp if  year==end_sample

twoway ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y>-50 & change_y<150  & sum_exp==1, ///
		msymbol(Oh) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mlw(medium) mcolor(gs5)) ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y>-50 & change_y<150  & sum_exp==2, ///
		msymbol(Oh) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mlw(thick) mcolor(gs5)) ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y>-50 & change_y<150  & sum_exp>=3, ///
		msymbol(O) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall)  mcolor(gs5)) ///
	(fpfit change_y dem_sum if e(sample) & year==end_sample, ///
		lcolor(gold*1) lw(.8) lpattern(solid)) ///
, xsize(11) ysize(6) scheme(s2mono)  graphregion(color(white) margin(small)) plotregion(margin(small)) ///
	xlabel(0(5)70, labsize(small)) yscale(range(-51 151)) ///
	ylabel(-50(50)150, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash) angle(0)) ///
	xtitle("Years Spent in Regime (Length of Treatment)", size(small)) ///
	ytitle("Total Change in Real GDP pc during Regime (in %)", size(small))  ///
	yline(0, lpattern(shortdash) lw(.4) lcolor(black)) xmtick(##5) ///
	text(125 0 "{bf:Freedom of Expression} (mean threshold)", size(medsmall) place(e)) ///
	legend(order(- "Number of Regime Changes" 1 2 3 4) label(1 "1") label(2 "2")  label(3 "3-4") ///
		label(4 "Fractional Polynomial") rows(1) size(small) ) 
graph export "$texfolder/Drilling_scatter_exp.eps", as(eps) replace

drop dem_sum end_sample yT start_sample y0 yT_ y0_ change_y change__y changey


*** FREEDOM OF ASSOCIATION
xtlogit dem_assoc_0 y dlpop ex if sample==1, fe
sort nwbcode year
by nwbcode: egen dem_sum=sum(sample) if e(sample) & dem_assoc_0==1
by nwbcode: egen end_sample=max(year) if e(sample) & dem_assoc_0==1
by nwbcode: egen start_sample=min(year) if e(sample) & dem_assoc_0==1
gen yT = ln(gdppc) if year==end_sample
by nwbcode: egen yT_ = mean(yT)
gen y0 = ln(gdppc) if year==start_sample
by nwbcode: egen y0_ = mean(y0)
gen changey = 100*(d.lngdppc) if e(sample) & dem_assoc_0==0
by nwbcode: egen change__y = mean(changey)
gen change_y = ((100*(yT_ - y0_)/dem_sum)/(change__y))*dem_sum

sort nwbcode year
by nwbcode: egen sum_assoc = sum(demevent_assoc) if sample==1
tab sum_assoc if  year==end_sample

twoway ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y<150 & change_y>-100 & sum_assoc==1, ///
		msymbol(Oh) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mlw(medium) mcolor(gs5)) ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y<150 & change_y>-100 & sum_assoc==2, ///
		msymbol(Oh) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mlw(thick) mcolor(gs5)) ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y<150 & change_y>-100 & sum_assoc>=3, ///
		msymbol(O) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mcolor(gs5)) ///
	(fpfit change_y dem_sum if e(sample) & year==end_sample, ///
		lcolor(purple*1.25) lw(.8) lpattern(solid)) ///
, xsize(11) ysize(6) scheme(s2mono)  graphregion(color(white) margin(small)) plotregion(margin(small)) ///
	xlabel(0(5)70, labsize(small))  ///
	ylabel(-100(50)150, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash) angle(0)) ///
	xtitle("Years Spent in Regime (Length of Treatment)", size(small)) ///
	ytitle("Total Change in Real GDP pc during Regime (in %)", size(small))  ///
	yline(0, lpattern(shortdash) lw(.4) lcolor(black)) xmtick(##5) ///
	text(125 0 "{bf:Freedom of Association} (mean threshold)", size(medsmall) place(e)) ///
	legend(order(- "Number of Regime Changes" 1 2 3 4) label(1 "1") label(2 "2")  label(3 "3-4") ///
		label(4 "Fractional Polynomial") rows(1) size(small) ) 
graph export "$texfolder/Drilling_scatter_assoc.eps", as(eps) replace

drop dem_sum end_sample yT start_sample y0 yT_ y0_ change_y change__y changey


*** FAIR ELECTIONS
xtlogit dem_fairel_0 y dlpop ex if sample==1, fe
sort nwbcode year
by nwbcode: egen dem_sum=sum(sample) if e(sample) & dem_fairel_0==1
by nwbcode: egen end_sample=max(year) if e(sample)
by nwbcode: egen start_sample=min(year) if e(sample)
gen yT = ln(gdppc) if year==end_sample
by nwbcode: egen yT_ = mean(yT)
gen y0 = ln(gdppc) if year==start_sample
by nwbcode: egen y0_ = mean(y0)
gen changey = 100*(d.lngdppc) if e(sample) & dem_assoc_0==0
by nwbcode: egen change__y = mean(changey)
gen change_y = ((100*(yT_ - y0_)/dem_sum)/(change__y))*dem_sum

sort nwbcode year
by nwbcode: egen sum_fairel = sum(demevent_fairel) if sample==1
tab sum_fairel if  year==end_sample

twoway ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y<150  & change_y>-50 & sum_fairel==1, ///
		msymbol(Oh) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mlw(medium) mcolor(gs5)) ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y<150  & change_y>-50 & sum_fairel==2, ///
		msymbol(Oh) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mlw(thick) mcolor(gs5)) ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y<150  & change_y>-50 & sum_fairel>=3, ///
		msymbol(O) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mcolor(gs5)) ///
	(fpfit change_y dem_sum if e(sample) & year==end_sample, ///
		lcolor(orange*1) lw(.8) lpattern(solid)) ///
, xsize(11) ysize(6) scheme(s2mono)  graphregion(color(white) margin(small)) plotregion(margin(small)) ///
	xlabel(0(5)70, labsize(small)) yscale(range(-51 151)) ///
	ylabel(-50(50)150, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash) angle(0)) ///
	xtitle("Years Spent in Regime (Length of Treatment)", size(small)) ///
	ytitle("Total Change in Real GDP pc during Regime (in %)", size(small))  ///
	yline(0, lpattern(shortdash) lw(.4) lcolor(black)) xmtick(##5) ///
	text(125 0 "{bf:Clean Elections} (mean threshold)", size(medsmall) place(e)) ///
	legend(order(- "Number of Regime Changes" 1 2 3 4) label(1 "1") label(2 "2")  label(3 "3-4") ///
		label(4 "Fractional Polynomial") rows(1) size(small) ) 
graph export "$texfolder/Drilling_scatter_fairel.eps", as(eps) replace

drop dem_sum end_sample yT start_sample y0 yT_ y0_ change_y change__y changey



*** RULE OF LAW
xtlogit dem_rolaw_0 y dlpop ex if sample==1, fe
sort nwbcode year
by nwbcode: egen dem_sum=sum(sample) if e(sample) & dem_rolaw_0==1
by nwbcode: egen end_sample=max(year) if e(sample)
by nwbcode: egen start_sample=min(year) if e(sample)
gen yT = ln(gdppc) if year==end_sample
by nwbcode: egen yT_ = mean(yT)
gen y0 = ln(gdppc) if year==start_sample
by nwbcode: egen y0_ = mean(y0)
gen changey = 100*(d.lngdppc) if e(sample) & dem_assoc_0==0
by nwbcode: egen change__y = mean(changey)
gen change_y = ((100*(yT_ - y0_)/dem_sum)/(change__y))*dem_sum

sort nwbcode year
by nwbcode: egen sum_rolaw = sum(demevent_rolaw) if sample==1
tab sum_rolaw if  year==end_sample

twoway ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y>-50  & change_y<150 & sum_rolaw==1, ///
		msymbol(Oh) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mlw(medium) mcolor(gs5)) ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y>-50  & change_y<150 & sum_rolaw==2, ///
		msymbol(Oh) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mlw(thick) mcolor(gs5)) ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y>-50  & change_y<150 & sum_rolaw>=3, ///
		msymbol(O) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mcolor(gs5)) ///
	(fpfit change_y dem_sum if e(sample) & year==end_sample, ///
		lcolor(midblue*1) lw(.8) lpattern(solid)) ///
, xsize(11) ysize(6) scheme(s2mono)  graphregion(color(white) margin(small)) plotregion(margin(small)) ///
	xlabel(0(5)70, labsize(small)) yscale(range(-51 151)) ///
	ylabel(-50(50)150, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash) angle(0)) ///
	xtitle("Years Spent in Regime (Length of Treatment)", size(small)) ///
	ytitle("Total Change in Real GDP pc during Regime (in %)", size(small))  ///
	yline(0, lpattern(shortdash) lw(.4) lcolor(black)) xmtick(##5) ///
	text(125 0 "{bf:Rule of Law} (mean threshold)", size(medsmall) place(e)) ///
	legend(order(- "Number of Regime Changes" 1 2 3 4) label(1 "1") label(2 "2")  label(3 "3-4") ///
		label(4 "Fractional Polynomial") rows(1) size(small) ) 
graph export "$texfolder/Drilling_scatter_rolaw.eps", as(eps) replace

drop dem_sum end_sample yT start_sample y0 yT_ y0_ change_y change__y changey

*** JUDICIAL CONSTRAINTS
xtlogit dem_jucon_0 y dlpop ex if sample==1, fe
sort nwbcode year
by nwbcode: egen dem_sum=sum(sample) if e(sample) & dem_jucon_0==1
by nwbcode: egen end_sample=max(year) if e(sample)
by nwbcode: egen start_sample=min(year) if e(sample)
gen yT = ln(gdppc) if year==end_sample
by nwbcode: egen yT_ = mean(yT)
gen y0 = ln(gdppc) if year==start_sample
by nwbcode: egen y0_ = mean(y0)
gen changey = 100*(d.lngdppc) if e(sample) & dem_assoc_0==0
by nwbcode: egen change__y = mean(changey)
gen change_y = ((100*(yT_ - y0_)/dem_sum)/(change__y))*dem_sum

sort nwbcode year
by nwbcode: egen sum_jucon = sum(demevent_jucon) if sample==1
tab sum_jucon if  year==end_sample


twoway ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y>-50  & change_y<150 & sum_jucon==1, ///
		msymbol(Oh) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mlw(medium) mcolor(gs5)) ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y>-50  & change_y<150 & sum_jucon==2, ///
		msymbol(Oh) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mlw(thick) mcolor(gs5)) ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y>-50  & change_y<150 & sum_jucon>=3, ///
		msymbol(O) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mcolor(gs5)) ///
	(fpfit change_y dem_sum if e(sample) & year==end_sample, ///
		lcolor(pink*1.25) lw(.8) lpattern(solid)) ///
, xsize(11) ysize(6) scheme(s2mono)  graphregion(color(white) margin(small)) plotregion(margin(small)) ///
	xlabel(0(5)70, labsize(small)) yscale(range(-51 151)) ///
	ylabel(-50(50)150, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash) angle(0)) ///
	xtitle("Years Spent in Regime (Length of Treatment)", size(small)) ///
	ytitle("Total Change in Real GDP pc during Regime (in %)", size(small))  ///
	yline(0, lpattern(shortdash) lw(.4) lcolor(black)) xmtick(##5) ///
	text(125 0 "{bf:Judicial Constraints on the Executive} (mean threshold)", size(medsmall) place(e)) ///
	legend(order(- "Number of Regime Changes" 1 2 3 4) label(1 "1") label(2 "2")  label(3 "3-4") ///
		label(4 "Fractional Polynomial") rows(1) size(small) ) 
graph export "$texfolder/Drilling_scatter_jucon.eps", as(eps) replace

drop dem_sum end_sample yT start_sample y0 yT_ y0_ change_y change__y changey


*** LEGISLATIVE CONSTRAINTS
xtlogit dem_legcon_0 y dlpop ex if sample==1, fe
sort nwbcode year
by nwbcode: egen dem_sum=sum(sample) if e(sample) & dem_legcon_0==1
by nwbcode: egen end_sample=max(year) if e(sample)
by nwbcode: egen start_sample=min(year) if e(sample)
gen yT = ln(gdppc) if year==end_sample
by nwbcode: egen yT_ = mean(yT)
gen y0 = ln(gdppc) if year==start_sample
by nwbcode: egen y0_ = mean(y0)
gen changey = 100*(d.lngdppc) if e(sample) & dem_assoc_0==0
by nwbcode: egen change__y = mean(changey)
gen change_y = ((100*(yT_ - y0_)/dem_sum)/(change__y))*dem_sum

sort nwbcode year
by nwbcode: egen sum_legcon = sum(demevent_legcon) if sample==1
tab sum_legcon if  year==end_sample

twoway ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y>-50 & change_y<150 & sum_legcon==1, ///
		msymbol(Oh) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mlw(medium) mcolor(gs5)) ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y>-50 & change_y<150 & sum_legcon==2, ///
		msymbol(Oh) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mlw(thick) mcolor(gs5)) ///
	(scatter change_y dem_sum if e(sample) & year==end_sample & change_y>-50 & change_y<150 & sum_legcon>=3, ///
		msymbol(O) msize(medlarge) mlab(wbcode) mlabposition(3) ///
		mlabsize(vsmall) mcolor(gs5)) ///
	(fpfit change_y dem_sum if e(sample) & year==end_sample, ///
		lcolor(emerald*1) lw(.8) lpattern(solid)) ///
, xsize(11) ysize(6) scheme(s2mono)  graphregion(color(white) margin(small)) plotregion(margin(small)) ///
	xlabel(0(5)70, labsize(small)) yscale(range(-51 151)) ///
	ylabel(-50(50)150, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash) angle(0)) ///
	xtitle("Years Spent in Regime (Length of Treatment)", size(small)) ///
	ytitle("Total Change in Real GDP pc during Regime (in %)", size(small))  ///
	yline(0, lpattern(shortdash) lw(.4) lcolor(black)) xmtick(##5) ///
	text(125 0 "{bf:Judicial Constraints on the Executive} (mean threshold)", size(medsmall) place(e)) ///
	legend(order(- "Number of Regime Changes" 1 2 3 4) label(1 "1") label(2 "2")  label(3 "3-4") ///
		label(4 "Fractional Polynomial") rows(1) size(small) ) 
graph export "$texfolder/Drilling_scatter_legcon.eps", as(eps) replace

drop dem_sum end_sample yT start_sample y0 yT_ y0_ change_y change__y changey

exit 


***
*** Sample makeup (Appendix)
***
tabstat year if sample==1, by(wbcode) statistics(min max N) columns(statistics)

gen dem_event_ld = (dem_ld==1 & l.dem_ld==0) if l.sample==1 & sample==1 & !missing(l.dem_ld)
gen rev_event_ld = (dem_ld==0 & l.dem_ld==1) if l.sample==1 & sample==1 & !missing(l.dem_ld)
gen dem_event_poly = (dem_poly==1 & l.dem_poly==0) if l.sample==1 & sample==1 & !missing(l.dem_poly)
gen rev_event_poly = (dem_poly==0 & l.dem_poly==1) if l.sample==1 & sample==1 & !missing(l.dem_poly)
gen dem_event_lib = (dem_lib==1 & l.dem_lib==0) if l.sample==1 & sample==1 & !missing(l.dem_lib)
gen rev_event_lib = (dem_lib==0 & l.dem_lib==1) if l.sample==1 & sample==1 & !missing(l.dem_lib)

preserve 
drop if sample!=1
sort wbcode year
collapse (firstnm) base_gdp = gdppc base_ld = v2x_libdem ///
	(lastnm) end_gdp = gdppc end_ld = v2x_libdem ///
	(sd) sd_gdp = gdppc sd_ld = v2x_libdem ///
	(sum) dem_event_ld rev_event_ld dem_event_poly rev_event_poly ///
	dem_event_lib rev_event_lib, by(wbcode)
list wbcode base_gdp end_gdp sd_gdp base_ld end_ld sd_ld  , separator(0) noobs
list wbcode dem_event_ld rev_event_ld dem_event_poly rev_event_poly dem_event_lib rev_event_lib ///
, separator(0) noobs
restore 
********************************************************************************
*** INCLUDED IN THE DRAFT
********************************************************************************

***
*** 
***

xtlogit dem_ld y if sample==1, fe iterate(1)
gen ln_gdppc = ln(gdppc)
graph box gdppc if e(sample), over(dem_ld) yscale(log) noout ylabel(500 1000 2500 5000 10000 20000 40000) ///
	name(g1, replace)
graph box gdppc if e(sample), over(dem_poly) yscale(log) noout ylabel(500 1000 2500 5000 10000 20000 40000) ///
	name(g2, replace)
graph box gdppc if e(sample), over(dem_lib) yscale(log) noout ylabel(500 1000 2500 5000 10000 20000 40000) ///
	name(g3, replace)

graph combine g1 g2 g3, name(G1) 


graph box wage, over(south) name(g1, replace) l1title(south) 
graph box wage, over(c_city) name(g2, replace) l1title(c_city) 
graph combine g1 g2, name(G1) 

stack wage south wage c_city, into(wage whatever) clear
label define whatever 0 zero 1 one 
label val whatever whatever 
label define _stack 1 south 2 c_city 
label val _stack _stack 

graph box wage, over(_stack) by(whatever) name(G2)


***
*** Conclusion: erosion of democratic institutions
***

use "$data/madd_ert_vdem_2024_helper2.dta", clear
gen data = 1 if !missing(v2x_libdem) & sample==1
sort nwbcode year
by nwbcode: gen sumdata = sum(data) 
gen demevent_ld = 1 if ((dem_ld==1 & l.dem_ld==0)) & sample==1
xtlogit dem_ld y if sample==1, fe iterate(1)
* Treat are those which experienced crossing the threshold: 66 countries
gen treat = (e(sample))

local z "_libdem _liberal _polyarchy _freexp_altinf _frassoc_thick el_frefair cl_rol _jucon lg_legcon _suffr _elecoff"
foreach k of local z{
	gen dv2x`k' = d.v2x`k'
}


gen dem_exp = (st_freexp_altinf>0)
gen dem_assoc = (st_frassoc_thick>0)
gen dem_fairel = (stel_frefair>0)

gen dem_rolaw = (stcl_rol>0)
gen dem_jucon = (st_jucon>0)
gen dem_legcon = (stlg_legcon>0)

tabstat year if treat==1, by(wbcode) statistics(N min max) columns(statistics)
* All treated countries have time series ending in 2018

* Suffrage and elected officials
tab v2x_suffr if sample==1
* 92% have university suffrage
gen dem_suffr= (v2x_suffr ==1)
replace dem_suffr = . if sample!=1
tab wbcode dem_suffr if sample==1
* Three countries

sum v2x_elecoff if sample==1
di .7720799-.25*.4064016
di .7720799+.25*.4064016

gen dem_elecoff= (v2x_elecoff > .7720799)
replace dem_elecoff = . if sample!=1
tab wbcode dem_elecoff if sample==1
* Twelve countries

drop dem_elecoff 
gen dem_elecoff= (v2x_elecoff > .8736803)
replace dem_elecoff = . if sample!=1
tab wbcode dem_elecoff if sample==1




sum dv2x* if treat==1 & year>=2009, separator(0)

tabstat dv2x_libdem if treat==1 & year>=2009, by(wbcode) statistics(sum) columns(statistics)
tabstat v2x_libdem if treat==1 & year==2009, by(wbcode) statistics(mean) columns(statistics)
tabstat dem_ld if treat==1, by(wbcode) statistics(sum) columns(statistics)
tabstat dem_ld if treat==1 & year==2018, by(wbcode) statistics(sum) columns(statistics)

* Polyarchy elements
tabstat dv2x_freexp_altinf if treat==1 & year>=2009, by(wbcode) statistics(N sum) columns(statistics)
tabstat v2x_freexp_altinf if treat==1 & year==2009, by(wbcode) statistics(mean) columns(statistics)
tabstat dem_exp if treat==1, by(wbcode) statistics(sum) columns(statistics)
tabstat dem_exp if treat==1 & year==2018, by(wbcode) statistics(sum) columns(statistics)

tabstat dv2x_frassoc_thick if treat==1 & year>=2009, by(wbcode) statistics(N sum) columns(statistics)
tabstat v2x_frassoc_thick if treat==1 & year==2009, by(wbcode) statistics(mean) columns(statistics)
tabstat dem_assoc if treat==1, by(wbcode) statistics(sum) columns(statistics)
tabstat dem_assoc if treat==1 & year==2018, by(wbcode) statistics(sum) columns(statistics)

tabstat dv2xel_frefair if treat==1 & year>=2009, by(wbcode) statistics(N sum) columns(statistics)
tabstat v2xel_frefair if treat==1 & year==2009, by(wbcode) statistics(mean) columns(statistics)
tabstat dem_fairel if treat==1, by(wbcode) statistics(sum) columns(statistics)
tabstat dem_fairel if treat==1 & year==2018, by(wbcode) statistics(sum) columns(statistics)


* Liberal component elements
tabstat dv2xcl_rol if treat==1 & year>=2009, by(wbcode) statistics(N sum) columns(statistics)
tabstat v2xcl_rol if treat==1 & year==2009, by(wbcode) statistics(mean) columns(statistics)
tabstat dem_rolaw if treat==1, by(wbcode) statistics(sum) columns(statistics)
tabstat dem_rolaw if treat==1 & year==2018, by(wbcode) statistics(sum) columns(statistics)

tabstat dv2x_jucon if treat==1 & year>=2009, by(wbcode) statistics(N sum) columns(statistics)
tabstat v2x_jucon if treat==1 & year==2009, by(wbcode) statistics(mean) columns(statistics)
tabstat dem_jucon if treat==1, by(wbcode) statistics(sum) columns(statistics)
tabstat dem_jucon if treat==1 & year==2018, by(wbcode) statistics(sum) columns(statistics)

tabstat dv2xlg_legcon if treat==1 & year>=2009, by(wbcode) statistics(N sum) columns(statistics)
tabstat v2xlg_legcon if treat==1 & year==2009, by(wbcode) statistics(mean) columns(statistics)
tabstat dem_legcon if treat==1, by(wbcode) statistics(sum) columns(statistics)
tabstat dem_legcon if treat==1 & year==2018, by(wbcode) statistics(sum) columns(statistics)

import excel "/Users/lezme/Dropbox/Facets of Democracy/Output/2021-12-02 Decadal decline analysis - conclusion.xlsx", sheet("Summary2") firstrow clear
sort order 

graph hbar (asis) declining declining_reg,  ///
	over(indicator, sort(order) relabel(6 "{bf:Liberal Democracy}")) ///
	ylabel(0(.1).7, labsize(small) glcolor(gs8) grid glwidth(.2) glp(shortdash) angle(0)) ///
	ytitle("Share of Treated Sample Countries", size(small)) yscale(range(-.01 .71)) ///
	scheme(s2mono)  graphregion(color(white) margin(small)) plotregion(margin(small)) ///
	bar(1, color(pink*1.5)) bar(2, color(midblue*.5)) ymtick(##5) ///
	legend(order(- "Share of Treated Countries Declining 2009-18" 1 2) rows(1) symxsize(*.35) ///
		label(1 "All") label(2 "In Regime") size(small)) ///
	text(.49 4 "-9.1%", place(w) just(right)) text(.54 9.5 "-12.8%", place(w) just(right) color(white)) ///
	text(.49 18.8 "-6.4%", place(w) just(right)) text(.54 24.3 "-7.7%", place(w) just(right) color(white)) ///
	text(.61 33.4 "-3.4%", place(w) just(right)) text(.58 38.8 "-4.3%", place(w) just(right) color(white)) ///
	text(.34 47.5 "-6.8%", place(w) just(right)) text(.42 53.3 "-16.4%", place(w) just(right) color(white)) ///
	text(.72 62 "-3.5%", place(w) just(right)) text(.63 67.3 "-4.0%", place(w) just(right) color(white)) ///
	text(.67 76.5 "-8.4%", place(w) just(right)) text(.6 82.4 "-8.7%", place(w) just(right) color(white)) ///
	text(.49 91.1 "-5.6%", place(w) just(right)) text(.57 96.2 "-13.3%", place(w) just(right) color(white))
graph export "$texfolder/decadal_decline.eps", as(eps) replace
	
	
